# This script is meant to be run from a facfAnalysisSetX directoy. It does the following:
# - integrate the force profile to get the free energy
# - output the file "GD.dat", containing free energies (G) and diffusion coefficients (D). This file is
#   later used in by the rp program to compute the permeability
# - collects the average constraint forces for each z-position in the single file "force.dat"

# usage:
# $ python computeG_collectGD.py 

import sys,string,math,os

o = open('force.dat','w') 

# read in force and diffusion values
f31 = open('31/D','r')
firstLine = f31.readline() # read first line
secondLine = f31.readline() # read second line
tokensFirst = string.split( firstLine )
tokensSecond = string.split( secondLine )
force31 = tokensFirst[0]
diff31 = tokensSecond[0]
f31.close()

f29 = open('29/D','r')
firstLine = f29.readline()
secondLine = f29.readline()
tokensFirst = string.split( firstLine )
tokensSecond = string.split( secondLine )
force29 = tokensFirst[0]
diff29 = tokensSecond[0]
f29.close()

f27 = open('27/D','r')
firstLine = f27.readline()
secondLine = f27.readline()
tokensFirst = string.split( firstLine )
tokensSecond = string.split( secondLine )
force27 = tokensFirst[0]
diff27 = tokensSecond[0]
f27.close()

f25 = open('25/D','r')
firstLine = f25.readline()
secondLine = f25.readline()
tokensFirst = string.split( firstLine )
tokensSecond = string.split( secondLine )
force25 = tokensFirst[0]
diff25 = tokensSecond[0]
f25.close()

f23 = open('23/D','r')
firstLine = f23.readline()
secondLine = f23.readline()
tokensFirst = string.split( firstLine )
tokensSecond = string.split( secondLine )
force23 = tokensFirst[0]
diff23 = tokensSecond[0]
f23.close()

f21 = open('21/D','r')
firstLine = f21.readline()
secondLine = f21.readline()
tokensFirst = string.split( firstLine )
tokensSecond = string.split( secondLine )
force21 = tokensFirst[0]
diff21 = tokensSecond[0]
f21.close()

f19 = open('19/D','r')
firstLine = f19.readline()
secondLine = f19.readline()
tokensFirst = string.split( firstLine )
tokensSecond = string.split( secondLine )
force19 = tokensFirst[0]
diff19 = tokensSecond[0]
f19.close()

f17 = open('17/D','r')
firstLine = f17.readline()
secondLine = f17.readline()
tokensFirst = string.split( firstLine )
tokensSecond = string.split( secondLine )
force17 = tokensFirst[0]
diff17 = tokensSecond[0]
f17.close()

f15 = open('15/D','r')
firstLine = f15.readline()
secondLine = f15.readline()
tokensFirst = string.split( firstLine )
tokensSecond = string.split( secondLine )
force15 = tokensFirst[0]
diff15 = tokensSecond[0]
f15.close()

f13 = open('13/D','r')
firstLine = f13.readline()
secondLine = f13.readline()
tokensFirst = string.split( firstLine )
tokensSecond = string.split( secondLine )
force13 = tokensFirst[0]
diff13 = tokensSecond[0]
f13.close()

f11 = open('11/D','r')
firstLine = f11.readline()
secondLine = f11.readline()
tokensFirst = string.split( firstLine )
tokensSecond = string.split( secondLine )
force11 = tokensFirst[0]
diff11 = tokensSecond[0]
f11.close()

f09 = open('09/D','r')
firstLine = f09.readline()
secondLine = f09.readline()
tokensFirst = string.split( firstLine )
tokensSecond = string.split( secondLine )
force09 = tokensFirst[0]
diff09 = tokensSecond[0]
f09.close()

f07 = open('07/D','r')
firstLine = f07.readline()
secondLine = f07.readline()
tokensFirst = string.split( firstLine )
tokensSecond = string.split( secondLine )
force07 = tokensFirst[0]
diff07 = tokensSecond[0]
f07.close()

f05 = open('05/D','r')
firstLine = f05.readline()
secondLine = f05.readline()
tokensFirst = string.split( firstLine )
tokensSecond = string.split( secondLine )
force05 = tokensFirst[0]
diff05 = tokensSecond[0]
f05.close()

f03 = open('03/D','r')
firstLine = f03.readline()
secondLine = f03.readline()
tokensFirst = string.split( firstLine )
tokensSecond = string.split( secondLine )
force03 = tokensFirst[0]
diff03 = tokensSecond[0]
f03.close()

f01 = open('01/D','r')
firstLine = f01.readline()
secondLine = f01.readline()
tokensFirst = string.split( firstLine )
tokensSecond = string.split( secondLine )
force01 = tokensFirst[0]
diff01 = tokensSecond[0]
f01.close()

# integrate force to get free energy differences
# and generate GD.dat file

dz = 0.2; # [nm] integration step

g31 = 0;
g29 = g31 + (float(force31)+float(force29))/2.0*dz;
g27 = g29 + (float(force29)+float(force27))/2.0*dz;
g25 = g27 + (float(force27)+float(force25))/2.0*dz;
g23 = g25 + (float(force25)+float(force23))/2.0*dz;
g21 = g23 + (float(force23)+float(force21))/2.0*dz;
g19 = g21 + (float(force21)+float(force19))/2.0*dz;
g17 = g19 + (float(force19)+float(force17))/2.0*dz;
g15 = g17 + (float(force17)+float(force15))/2.0*dz;
g13 = g15 + (float(force15)+float(force13))/2.0*dz;
g11 = g13 + (float(force13)+float(force11))/2.0*dz;
g09 = g11 + (float(force11)+float(force09))/2.0*dz;
g07 = g09 + (float(force09)+float(force07))/2.0*dz;
g05 = g07 + (float(force07)+float(force05))/2.0*dz;
g03 = g05 + (float(force05)+float(force03))/2.0*dz;
g01 = g03 + (float(force03)+float(force01))/2.0*dz;

gd = open('GD.dat','w')
print >> gd, g31, "\t", diff31
print >> gd, g29, "\t", diff29
print >> gd, g27, "\t", diff27
print >> gd, g25, "\t", diff25
print >> gd, g23, "\t", diff23
print >> gd, g21, "\t", diff21
print >> gd, g19, "\t", diff19
print >> gd, g17, "\t", diff17
print >> gd, g15, "\t", diff15
print >> gd, g13, "\t", diff13
print >> gd, g11, "\t", diff11
print >> gd, g09, "\t", diff09
print >> gd, g07, "\t", diff07
print >> gd, g05, "\t", diff05
print >> gd, g03, "\t", diff03
print >> gd, g01, "\t", diff01



# print out "force.dat" file
# this is not strictly necessary  
print >> o, "-3.1\t", force31
print >> o, "-2.9\t", force29
print >> o, "-2.7\t", force27 
print >> o, "-2.5\t", force25 
print >> o, "-2.3\t", force23 
print >> o, "-2.1\t", force21 
print >> o, "-1.9\t", force19 
print >> o, "-1.7\t", force17 
print >> o, "-1.5\t", force15 
print >> o, "-1.3\t", force13 
print >> o, "-1.1\t", force11 
print >> o, "-0.9\t", force09 
print >> o, "-0.7\t", force07 
print >> o, "-0.5\t", force05 
print >> o, "-0.3\t", force03 
print >> o, "-0.1\t", force01

o.close()
